Error Correction in Ancestry Classification

ABSTRACT

Error correction in ancestry classification includes obtaining, from a classifier, initial ancestry classifications associated with portions of two phased haplotypes of a chromosome pair of an individual; performing error correction on an initial ancestry classification, including detecting a phasing error in the initial ancestry classifications; and outputting a corrected ancestry classification in which the phasing error is corrected.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.18/157,595, filed on Jan. 20, 2023; which is a continuation of U.S.patent application Ser. No. 16/844,758, filed on Apr. 9, 2020; which isa continuation of U.S. patent application Ser. No. 14/938,111, filed onNov. 11, 2015 and issued as U.S. Pat. No. 10,658,071 on May 19, 2020;which is a continuation of U.S. patent application Ser. No. 13/801,056,filed on Mar. 13, 2013 and issued as U.S. Pat. No. 9,213,947 on Dec. 15,2015; which claims the benefit of U.S. Provisional Patent ApplicationNo. 61/724,236, filed on Nov. 8, 2012, and U.S. Provisional PatentApplication No. 61/724,228, filed on Nov. 8, 2012. The foregoingapplications are incorporated herein by reference in their entireties.

SEQUENCE LISTING STATEMENT

A computer readable form of the Sequence Listing is filed with thisapplication by electronic submission and is incorporated into thisapplication by reference in its entirety. The Sequence Listing iscontained in the file created on Jul. 14, 2023, having the file name“22-1266-US-CON5_Sequence-Listing.xml” and is 4,096 bytes in size.

BACKGROUND OF THE INVENTION

Ancestry deconvolution refers to identifying the ancestral origin ofchromosomal segments in individuals. Ancestry deconvolution in admixedindividuals (i.e., individuals whose ancestors such as grandparents arefrom different regions) is straightforward when the ancestralpopulations considered are sufficiently distinct (e.g., one grandparentis from Europe and another from Asia). To date, however, existingapproaches are typically ineffective at distinguishing between closelyrelated populations (e.g., within Europe). Moreover, due to theircomputational complexity, most existing methods for ancestrydeconvolution are unsuitable for application in large-scale settings,where the reference panels used contain thousands of individuals.

BRIEF DESCRIPTION OF THE DRAWINGS

Various embodiments of the invention are disclosed in the followingdetailed description and the accompanying drawings.

FIG. 1 is a functional diagram illustrating a programmed computer systemfor performing the pipelined ancestry prediction process in accordancewith some embodiments.

FIG. 2 is a block diagram illustrating an embodiment of an ancestryprediction platform.

FIG. 3 is an architecture diagram illustrating an embodiment of anancestry prediction system.

FIG. 4 is a flowchart illustrating an embodiment of a process forancestry prediction.

FIG. 5A illustrates an example of a section of unphased genotype data.

FIG. 5B illustrates an example of two sets of phased genotype data.

FIG. 6 is a flowchart illustrating an embodiment of a process forperforming out-of-sample phasing.

FIG. 7 is a diagram illustrating an example of a predetermined referencehaplotype graph that is built based on a reference collection ofgenotype data.

FIGS. 8A-8B are diagrams illustrating embodiments of modified haplotypegraphs.

FIG. 9 is a diagram illustrating an embodiment of a compressed haplotypegraph with segments.

FIG. 10 is a diagram illustrating an embodiment of a dynamic Bayesiannetwork used to implement trio-based phasing.

FIG. 11 is a flowchart illustrating an embodiment of a process toperform trio-based phasing.

FIG. 12 is a flowchart illustrating an embodiment of a localclassification process.

FIG. 13 is a diagram illustrating how a set of reference data points isclassified into two classes by a binary SVM.

FIG. 14 is a diagram illustrating possible errors in an exampleclassification result.

FIG. 15 is a graph illustrating embodiments of Hidden Markov Models usedto model phasing errors.

FIG. 16 shows two example graphs corresponding to basic HMMs used tomodel phasing errors of two example haplotypes of a chromosome.

FIG. 17 is an example graph of a Pair Hidden Markov Model (PHMM).

FIG. 18 is an example graph of an Autoregressive Pair Hidden MarkovModel (APHMM).

FIG. 19 is an example data table displaying the emission parameters.

FIGS. 20A-20D are example ancestry assignment plots illustratingdifferent results that are obtained using different techniques.

FIG. 21 is a table comparing the predictive accuracies of ancestryassignments with and without error correction.

FIG. 22A illustrates example reliability plots for East Asian populationand East European population before recalibration.

FIG. 22B illustrates example reliability plots for East Asian populationand East European population after recalibration.

FIG. 23A is a flowchart illustrating an embodiment of a label clusteringprocess.

FIG. 23B is an example illustrating process 2300 of FIG. 23A.

FIG. 24 is a flowchart illustrating an embodiment of a process fordisplaying ancestry information.

FIG. 25 is a diagram illustrating an embodiment of a regional view ofancestry composition information for an individual.

FIG. 26 is a diagram illustrating an embodiment of an expanded view ofancestry composition information for an individual.

FIG. 27 is a diagram illustrating an embodiment of a further expandedview of ancestry composition information for an individual.

FIG. 28 is a diagram illustrating an embodiment of an inheritance view.

FIGS. 29-30 are diagrams illustrating embodiments of achromosome-specific view.

DETAILED DESCRIPTION

The invention can be implemented in numerous ways, including as aprocess; an apparatus; a system; a composition of matter; a computerprogram product embodied on a computer readable storage medium; and/or aprocessor, such as a processor configured to execute instructions storedon and/or provided by a memory coupled to the processor. In thisspecification, these implementations, or any other form that theinvention may take, may be referred to as techniques. In general, theorder of the steps of disclosed processes may be altered within thescope of the invention. Unless stated otherwise, a component such as aprocessor or a memory described as being configured to perform a taskmay be implemented as a general component that is temporarily configuredto perform the task at a given time or a specific component that ismanufactured to perform the task. As used herein, the term ‘processor’refers to one or more devices, circuits, and/or processing coresconfigured to process data, such as computer program instructions.

A detailed description of one or more embodiments of the invention isprovided below along with accompanying figures that illustrate theprinciples of the invention. The invention is described in connectionwith such embodiments, but the invention is not limited to anyembodiment. The scope of the invention is limited only by the claims andthe invention encompasses numerous alternatives, modifications andequivalents. Numerous specific details are set forth in the followingdescription in order to provide a thorough understanding of theinvention. These details are provided for the purpose of example and theinvention may be practiced according to the claims without some or allof these specific details. For the purpose of clarity, technicalmaterial that is known in the technical fields related to the inventionhas not been described in detail so that the invention is notunnecessarily obscured.

A pipelined ancestry deconvolution process to predict an individual'sancestry based on genetic information is disclosed. Unphased genotypedata associated with the individual's chromosomes is received and phasedto generate phased haplotype data. In some embodiments, dynamicprogramming that does not require the unphased genotype data to beincluded in the reference data is implemented to facilitate phasing. Thephased data is divided into segments, which are classified as beingassociated with specific ancestries. The classification is performedusing a learning machine in some embodiments. The classification outputundergoes an error correction process to reduce noise and correct forany phasing errors (also referred to as switch errors) and/or correlatedclassification errors. The error corrected output is optionallyrecalibrated, and ancestry labels are optionally clustered according toa geographical hierarchy to be displayed to the user.

In some embodiments, genotype data comprising gene sequences and/orgenetic markers is used to represent an individual's genome. Examples ofsuch genetic markers include Single Nucleotide Polymorphisms (SNPs),which are points along the genome, each corresponding to two or morecommon variations; Short Tandem Repeats (STRs), which are repeatedpatterns of two or more repeated nucleotide sequences adjacent to eachother; and Copy-Number Variants (CNVs), which include longer sequencesof deoxyribonucleic acid (DNA) that could be present in varying numbersin different individuals. Although SNP-based genotype data is describedextensively below for purposes of illustration, the technique is alsoapplicable to other types of genotype data such as STRs and CNVs. Asused herein, a haplotype refers to DNA on a single chromosome of achromosome pair. Haplotype data representing a haplotype can beexpressed as a set of markers (e.g., SNPs, STRs, CNVs, etc.) or a fullDNA sequence set.

FIG. 1 is a functional diagram illustrating a programmed computer systemfor performing the pipelined ancestry prediction process in accordancewith some embodiments. Computer system 100, which includes varioussubsystems as described below, includes at least one microprocessorsubsystem (also referred to as a processor or a central processing unit(CPU)) 102. For example, processor 102 can be implemented by asingle-chip processor or by multiple processors. In some embodiments,processor 102 is a general purpose digital processor that controls theoperation of the computer system 100. Using instructions retrieved frommemory 110, the processor 102 controls the reception and manipulation ofinput data, and the output and display of data on output devices (e.g.,display 118). In some embodiments, processor 102 includes and/or is usedto provide phasing, local classification, error correction,recalibration, and/or label clustering as described below.

Processor 102 is coupled bi-directionally with memory 110, which caninclude a first primary storage, typically a random access memory (RAM),and a second primary storage area, typically a read-only memory (ROM).As is well known in the art, primary storage can be used as a generalstorage area and as scratch-pad memory, and can also be used to storeinput data and processed data. Primary storage can also storeprogramming instructions and data, in the form of data objects and textobjects, in addition to other data and instructions for processesoperating on processor 102. Also as is well known in the art, primarystorage typically includes basic operating instructions, program code,data, and objects used by the processor 102 to perform its functions(e.g., programmed instructions). For example, memory 110 can include anysuitable computer-readable storage media, described below, depending onwhether, for example, data access needs to be bi-directional oruni-directional. For example, processor 102 can also directly and veryrapidly retrieve and store frequently needed data in a cache memory (notshown).

A removable mass storage device 112 provides additional data storagecapacity for the computer system 100, and is coupled eitherbi-directionally (read/write) or uni-directionally (read only) toprocessor 102. For example, storage 112 can also includecomputer-readable media such as magnetic tape, flash memory, PC-CARDS,portable mass storage devices, holographic storage devices, and otherstorage devices. A fixed mass storage 120 can also, for example, provideadditional data storage capacity. The most common example of massstorage 120 is a hard disk drive. Mass storage 112, 120 generally storeadditional programming instructions, data, and the like that typicallyare not in active use by the processor 102. It will be appreciated thatthe information retained within mass storage 112 and 120 can beincorporated, if needed, in standard fashion as part of memory 110(e.g., RAM) as virtual memory.

In addition to providing processor 102 access to storage subsystems, bus114 can also be used to provide access to other subsystems and devices.As shown, these can include a display monitor 118, a network interface116, a keyboard 104, and a pointing device 106, as well as an auxiliaryinput/output device interface, a sound card, speakers, and othersubsystems as needed. For example, the pointing device 106 can be amouse, stylus, track ball, or tablet, and is useful for interacting witha graphical user interface.

The network interface 116 allows processor 102 to be coupled to anothercomputer, computer network, or telecommunications network using anetwork connection as shown. For example, through the network interface116, the processor 102 can receive information (e.g., data objects orprogram instructions) from another network or output information toanother network in the course of performing method/process steps.Information, often represented as a sequence of instructions to beexecuted on a processor, can be received from and outputted to anothernetwork. An interface card or similar device and appropriate softwareimplemented by (e.g., executed/performed on) processor 102 can be usedto connect the computer system 100 to an external network and transferdata according to standard protocols. For example, various processembodiments disclosed herein can be executed on processor 102, or can beperformed across a network such as the Internet, intranet networks, orlocal area networks, in conjunction with a remote processor that sharesa portion of the processing. Additional mass storage devices (not shown)can also be connected to processor 102 through network interface 116.

An auxiliary I/O device interface (not shown) can be used in conjunctionwith computer system 100. The auxiliary I/O device interface can includegeneral and customized interfaces that allow the processor 102 to sendand, more typically, receive data from other devices such asmicrophones, touch-sensitive displays, transducer card readers, tapereaders, voice or handwriting recognizers, biometrics readers, cameras,portable mass storage devices, and other computers.

In addition, various embodiments disclosed herein further relate tocomputer storage products with a computer readable medium that includesprogram code for performing various computer-implemented operations. Thecomputer-readable medium is any data storage device that can store datawhich can thereafter be read by a computer system. Examples ofcomputer-readable media include, but are not limited to, all the mediamentioned above: magnetic media such as hard disks, floppy disks, andmagnetic tape; optical media such as CD-ROM disks; magneto-optical mediasuch as optical disks; and specially configured hardware devices such asapplication-specific integrated circuits (ASICs), programmable logicdevices (PLDs), and ROM and RAM devices. Examples of program codeinclude both machine code, as produced, for example, by a compiler, orfiles containing higher level code (e.g., script) that can be executedusing an interpreter.

The computer system shown in FIG. 1 is but an example of a computersystem suitable for use with the various embodiments disclosed herein.Other computer systems suitable for such use can include additional orfewer subsystems. In addition, bus 114 is illustrative of anyinterconnection scheme serving to link the subsystems. Other computerarchitectures having different configurations of subsystems can also beutilized.

FIG. 2 is a block diagram illustrating an embodiment of an ancestryprediction platform. In this example, a user uses a client device 202 tocommunicate with an ancestry prediction system 206 via a network 204.Examples of device 202 include a laptop computer, a desktop computer, asmart phone, a mobile device, a tablet device or any other computingdevice. Ancestry prediction system 206 is used to perform a pipelinedprocess to predict ancestry based on a user's genotype information.Ancestry prediction system 206 can be implemented on a networkedplatform (e.g., a server or cloud-based platform, a peer-to-peerplatform, etc.) that supports various applications. For example,embodiments of the platform perform ancestry prediction and provideusers with access (e.g., via appropriate user interfaces) to theirpersonal genetic information (e.g., genetic sequence information and/orgenotype information obtained by assaying genetic materials such asblood or saliva samples) and predicted ancestry information. In someembodiments, the platform also allows users to connect with each otherand share information. Device 100 can be used to implement 202 or 206.

In some embodiments, DNA samples (e.g., saliva, blood, etc.) arecollected from genotyped individuals and analyzed using DNA microarrayor other appropriate techniques. The genotype information is obtained(e.g., from genotyping chips directly or from genotyping services thatprovide assayed results) and stored in database 208 and is used bysystem 206 to make ancestry predictions. Reference data, includinggenotype data of unadmixed individuals (e.g., individuals whoseancestors came from the same region), simulated data (e.g., results ofmachine-based processes that simulate biological processes such asrecombination of parents' DNA), pre-computed data (e.g., a precomputedreference haplotype graph used in out-of-sample phasing) and the likecan also be stored in database 208 or any other appropriate storageunit.

FIG. 3 is an architecture diagram illustrating an embodiment of anancestry prediction system. System 300 can be used to implement 206 ofFIG. 2 , and can be implemented using system 100 of FIG. 1 . Theprocessing pipeline of system 300 includes a phasing module 302, a localclassification module 304, and an error correction module 306. Thesemodules form a predictive engine that makes predictions about therespective ancestries that correspond to the individual's chromosomeportions. Optionally, a recalibration module 308 and/or a labelclustering module 310 can also be included to refine the output of thepredictive engine.

The input to phasing module 302 comprises unphased genotype data, andthe output of the phasing module comprises phased genotype data (e.g.,two sets of haplotype data). In some embodiments, phasing module 302performs out-of-sample phasing where the unphased genotype data beingphased is not included in the reference data used to perform phasing.The phased genotype data is input into local classification module 304,which outputs predicted ancestry information associated with the phasedgenotype data. In some embodiments, the phased genotype data issegmented, and the predicted ancestry information includes one or moreancestry predictions associated with the segments. The posteriorprobabilities associated with the predictions are also optionallyoutput. The predicted ancestry information is sent to error correctionmodule 306, which averages out noise in the predicted ancestryinformation and corrects for phasing errors introduced by the phasingmodule and/or correlated prediction errors introduced by the localclassification module. The output of the error correction module can bepresented to the user (e.g., via an appropriate user interface).Optionally, the error correction module sends its output (e.g., errorcorrected posterior probabilities) to a recalibration module 308, whichrecalibrates the output to establish confidence levels based on theerror corrected posterior probabilities. Also optionally, the calibratedconfidence levels are further sent to label clustering module 310 toidentify appropriate ancestry assignments that meet a confidence levelrequirement.

The modules described above can be implemented as software componentsexecuting on one or more processors, as hardware such as programmablelogic devices and/or Application Specific Integrated Circuits designedto perform certain functions or a combination thereof. In someembodiments, the modules can be embodied by a form of software productswhich can be stored in a nonvolatile storage medium (such as opticaldisk, flash storage device, mobile hard disk, etc.), including a numberof instructions for making a computer device (such as personalcomputers, servers, network equipment, etc.) implement the methodsdescribed in the embodiments of the present application. The modules maybe implemented on a single device or distributed across multipledevices. The functions of the modules may be merged into one another orfurther split into multiple sub-modules.

In addition to being a part of the pipelined ancestry predictionprocess, the modules and their outputs can be used in otherapplications. For example, the output of the phasing module can be usedto identify familial relatives of individuals in the reference database.

FIG. 4 is a flowchart illustrating an embodiment of a process forancestry prediction. Process 400 initiates at 402, when unphasedgenotype data associated with one or more chromosomes of an individualis obtained. The unphased genotype data can be received from a datasource such as a database or a genotyping service, or obtained by userupload. At 404, the unphased genotype data is phased using anout-of-sample technique to generate two sets of phased haplotype data.Each set of phased haplotype data corresponds to the DNAs the individualinherited from one biological parent. At 406, a learning machine (e.g.,a support vector machine (SVM)) is used to classify portions of the twosets of haplotype data as being associated with specific ancestriesrespectively, and generate ancestry classification results. At 408,errors in the results of the ancestry classification are corrected. Insome embodiments, error correction removes noise, corrects phasingerrors and/or correlated prediction errors. Optionally, at 410, theerror corrected predicted ancestry information is recalibrated toestablish confidence levels. Optionally, at 412, the recalibratedconfidence levels and their associated ancestry assignments areclustered as appropriate to identify ancestry assignments that meet aconfidence level requirement. Optionally, at 414, the resultingconfidence levels and their associated ancestry assignments are storedto a database and/or output to another application (e.g., an applicationthat analyzes the results and/or displays predicted ancestry informationto users).

Details of the modules and their operations are described below.

Phasing

At a given gene locus on a pair of autosomal chromosomes, a diploidorganism (e.g., a human being) inherits one allele of the gene from themother and another allele of the gene from the father. At a heterozygousgene locus, two parents contribute different alleles (e.g., one A andone C). Without additional processing, it is impossible to tell whichparent contributed which allele. Such genotype data that is notattributed to a particular parent is referred to as unphased genotypedata. Typically, initial genotype readings obtained from genotypingchips manufactured by companies such as Illumina® are in an unphasedform.

FIG. 5A illustrates an example of a section of unphased genotype data.Genotype data section 502 includes genotype calls at known SNP locationsof a chromosome pair. The process of phasing is to split a stretch ofunphased genotype calls such as 502 into two sets of phased genotypedata (also referred to as haplotype data) attributed to a particularparent. Phasing is needed for identifying ancestry from each parent andclassifying haplotypes from different ancestral origins. Further, aspecific marker alone tends not to offer good ancestral (e.g.,geographical or ethic) specificity, but a run of multiple markers canoffer better specificity. For example, a particular SNP of “A” is notvery informative with respect to the ancestry origin of the section ofDNA, but a haplotype of a longer stretch (e.g., “ACGA”) starting at aspecific location can be highly correlated with Northern Europeanancestry.

FIG. 5B illustrates an example of two sets of phased genotype data. Inthis example, phased genotype data (i.e., haplotype data) 504 and 506 isobtained from unphased genotype data 502 based on statisticaltechniques. Haplotype block 504 (“ACGT”) is determined to be attributedto (i.e., inherited from) one parent, and haplotype block 506 (“AACC”)is determined to be attributed to another parent.

Population-Based Phasing

Phasing is often done using statistical techniques. Such techniques arealso referred to as population-based phasing because genotype data froma reference collection of a population of individuals (e.g., a fewhundred to a thousand) is analyzed. BEAGLE is a commonly usedpopulation-based phasing technique. It makes statistical determinationsbased on the assumption that certain blocks of haplotypes are inheritedin blocks and therefore shared amongst individuals. For example, if thegenotype data of a sample population comprising many individuals shows acommon pattern of “?A ?C ?G ?T” (where “?” can be any other allele),then the block “ACGT” is likely to be a common block of haplotypes thatis present in these individuals. The population-based phasing techniquewould therefore identify the block “ACGT” as coming from one parentwhenever “?A ?C ?G ?T” is present in the genotype data. Because BEAGLErequires that the genotype data being analyzed be included in thereference collection, the technique is referred to as in-sample phasing.

In-sample phasing is often computationally inefficient. Phasing of alarge database of a user's genome (e.g., 100,000 or more) can take manydays, and it can take just as long whenever a new user has to be addedto the database since the technique would recompute the full set of data(including the new user's data). There can also be mistakes duringin-sample phasing. One type of mistake, referred to as phasing errors orswitch errors, occurs where a section of the chromosome is in factattributed to one parent but is misidentified as attributed to anotherparent. Switch errors can occur when a stretch of genotype data is notcommon in the reference population. For example, suppose that a parentactually contributed the haplotype of “ACCC” and another parent actuallycontributed the haplotype of “AAGT” to genotype 502. Because the block“ACGT” is common in the reference collection and “ACCC” has neverappeared in the reference collection, the technique attributes “ACGT”and “AACC” to two parents respectively, resulting in a switch error.

Embodiments of the phasing technique described below permitout-of-sample population-based phasing. In out-of-sample phasing, whengenotype data of a new individual needs to be phased, the genotype datais not necessarily immediately combined with the reference collection toobtain phasing for this individual. Instead, a precomputed datastructure such as a predetermined reference haplotype graph is used tofacilitate a dynamic programming based process that quickly phases thegenotype data. For example, given the haplotype graph and unphased data,the likely sequence of genotype data can be solved using the Viterbialgorithm. This way, on a platform with a large number of users forminga large reference collection (e.g., at least 100,000 individuals), whena new individual signs up with the service and provides his/her genotypedata, the platform is able to quickly phase the genotype data withouthaving to recompute the common haplotypes of the existing users plus thenew individual.

FIG. 6 is a flowchart illustrating an embodiment of a process forperforming out-of-sample phasing. Process 600 can be performed on asystem such as 100 or 206, and can be used to implement phasing module302.

At 602, unphased genotype data of the individual is obtained. In someembodiments, the unphased genotype data such as sequence data 502 isreceived from a database, a genotyping service, or as an upload by auser of a platform such as 100.

At 604, the unphased genotype data is processed using dynamicprogramming to determine phased data, i.e., sets of likely haplotypes.The processing requires a reference population and is therefore referredto as population-based phasing. In some embodiments, the dynamicprogramming relies on a predetermined reference haplotype graph. Thepredetermined haplotype graph is precomputed without referencing theunphased genotype data of the individual. Thus, the unphased genotypedata is said to be out-of-sample with respect to a collection ofreference genotype data used to compute the predetermined referencehaplotype graph. In other words, if the unphased genotype data is from anew user whose genotype data is not already included in the referencegenotype data and therefore is not incorporated into the predeterminedreference haplotype graph, it is not necessary to include the unphasedgenotype data from the new user in the reference genotype data andrecompute the reference haplotype graph. Details of dynamic programmingand the predetermined reference haplotype graph are described below.

At 606, trio-based phasing is optionally performed to improve upon theresults from population-based phasing. As used herein, trio-basedphasing refers to phasing by accounting for the genotyping data of oneor more biological parents of the individual.

At 608, the likely haplotype data is output to be stored to a databaseand/or processed further. In some embodiments, the likely haplotype datais further processed by a local classifier as shown in FIG. 3 forancestry prediction purposes.

The likely haplotype data can also be used in other applications, suchas being compared with haplotype data of other individuals in a databaseto identify the amount of DNA shared among individuals, therebydetermining people who are related to each other and/or people belongingto the same population groups.

In some embodiments, the dynamic programming process performed in step604 uses a predetermined reference haplotype graph to examine possiblesequences of haplotypes that could be combined to generate the unphasedgenotype data, and determine the most likely sequences of haplotypes.Given a collection of binary strings of length L, a haplotype graph is aprobabilistic deterministic finite automaton (DFA) defined over adirected acyclic graph. The nodes of the multigraph are organized intoL+1 levels (numbered from 0 to L), such that level 0 has a single noderepresenting the source (i.e., initial state) of the DFA and level L hasa single node representing the sink (i.e., accepting state) of the DFA.Every directed edge in the multigraph connects a node from some level ito a node in level (i+1) and is labeled with either 0 or 1. Every nodeis reachable from the source and has a directed path to the sink. Foreach path through the haplotype graph from the source to the sink, theconcatenation of the labels on the edges traversed by the path is abinary string of length L. Semantically, paths through the graphrepresent haplotypes over a genomic region comprising L biallelicmarkers (assuming an arbitrary binary encoding of the alleles at eachsite). A probability distribution over the set of haplotypes included ina haplotype graph can be defined by associating a conditionalprobability with each edge (such that the sum of the probabilities ofthe outgoing edges for each node is equal to 1), and generated bystarting from the initial state at level 0, and choosing successorstates by following random outgoing edges according to their assignedconditional probabilities.

FIG. 7 is a diagram illustrating an example of a predetermined referencehaplotype graph that is built based on a reference collection ofgenotype data (e.g., population-based data). In this example, thereference collection of genotype data includes a set of L geneticmarkers (e.g., SNPs). Haplotype graph 700 is a Directed Acyclic Graph(DAG) having nodes (e.g., 704) and edges (e.g., 706). The haplotypegraph starts with a single node (the “begin state”) and ends on a singlenode (the “accepting state”), and the intermediate nodes correspond tothe states of the markers at respective gene loci. There are a total ofL+1 levels of nodes from left to right. An edge, e, represents the setof haplotypes whose path from the initial node to the terminating nodeof the graph traverses e. The possible paths define the haplotype spaceof possible genotype sequences. For example, in haplotype graph 700, apossible path 702 corresponds to the genotype sequence “GTTCAC”. Thereare four possible paths/genotype sequences in the haplotype space shownin this diagram (“ACGCGC,” “ACTTAC,” “GTTCAC,” and “GTTTGG”).

Each edge is associated with a probability computed based on thereference collection of genotype data. In this example, a collection ofgenotype data is comprised of genotype data from 1000 individuals, ofwhich 400 have the “A” allele at the first locus, and 600 have the “G”allele at the first locus. Accordingly, the probability associated withedge 708 is 400/1000 and the probability associated with edge 710 is600/1000. All of the first 400 individuals have the “C” allele at thesecond locus, giving edge 712 a probability of 400/400. All of the next600 individuals who had the “G” allele at the first locus have the “T”allele at the second locus, giving edge 714 a probability of 600/600,and so on. The probabilities associated with the respective edges arelabeled in the diagram. The probability associated with a specific pathis expressed as the product of the probabilities associated with theedges included in the path. For example, the probability associated withpath 702 is computed as:

${P(h)} = {{\left( \frac{600}{1000} \right)\left( \frac{600}{600} \right)\left( \frac{600}{600} \right)\left( \frac{50}{600} \right)\left( \frac{350}{350} \right)\left( \frac{450}{450} \right)} = 0.05}$

The dynamic programming process searches the haplotype graph forpossible paths, selecting two paths h₁ and h₂ for which the product oftheir associated probabilities is maximized, subject to the constraintthat when the two paths are combined, the alleles at each locus mustmatch the corresponding alleles in the unphased genotype data (g). Thefollowing expression is used in some cases to characterize the process:

maximize P(h ₁)P(h ₂), subject to h ₁ +h ₂ =g

For out-of-sample phasing, the reference haplotype graph is built onceand reused to identify possible haplotype paths that correspond to theunphased genotype data of a new individual (a process also referred toas “threading” the new individual's haplotype along the graph). Theindividual's genotype data sometimes does not correspond to any existingpath in the graph (e.g., the individual has genotype sequences that areunique and not included in the reference population), and thereforecannot be successfully threaded based on existing paths of the referencehaplotype graph. To cope with the possibility of a non-existent path,several modifications are made to the reference haplotype graph tofacilitate the out-of-sample phasing process.

FIGS. 8A-8B are diagrams illustrating embodiments of modified haplotypegraph used for out-of-sample, population-based phasing. In theseexamples, modified reference haplotype graphs 800 and 850 are based ongraph 700. Unlike graph 700, which is based on exact readings ofgenotype sequences of the reference individuals, the modified graphspermit recombination and genotyping errors and include modifications(e.g., extra edges) that account for recombination and genotypingerrors.

Recombination is one reason to extend graph 700 for out-of-samplephasing. As used herein, recombination refers to the switching of ahaplotype along one path to a different path. Recombination can happenwhen segments of parental chromosomes cross over during meiosis. In someembodiments, reference haplotype graph 700 is extended to account forthe possibility of recombination/path switching. Recombination eventsare modeled by allowing a new haplotype state to be selected(independent of the previous haplotype state) with probability τ at eachlevel of the haplotype graph. By default, τ≈0.00448, which is anestimate of the probability of recombination between adjacent sites,assuming 500,000 uniformly spaced markers, a genome length of 37.5Morgans, and 30 generations since admixture. Referring to the example ofFIG. 8A, suppose the new individual's unphased genotype data is “AG, CT,TT, TT, GG, GG,” (SEQ ID NO: 1) which cannot be split into twohaplotypes by threading along existing paths in graph 700. The modifiedreference haplotype graph 800 permits recombination by includingadditional edges representing recombination (e.g., edge 804) so that newpaths can be formed along these edges. In this example, the unphasedgenotype data can map onto two paths corresponding to haplotypes“ACTTGG” and “GTTTGG”, the former being a new path due to recombinationwith a recombination occurring between “C” and “T” along edge 804. τ isassociated with edge 804 and used to compute the probability of the paththrough 804.

Genotyping error is another reason to extend graph 700 for out-of-samplephasing. Genotyping errors can occur because the genotyping technologyis imperfect and can make false readings. The rate of genotyping errorfor a given technology (e.g., a particular genotyping chip) can beobtained from the manufacturer. In some embodiments, when the search forpossible paths for a new individual cannot be done according to theexisting reference graph, the existing reference haplotype graph isextended to account for the possibility of genotyping errors. Forexample, suppose the new individual's unphased genotype data is “AG, CT,GG, CT, GG, CG,” (SEQ ID NO: 2) which cannot be split into twohaplotypes by threading along existing paths in graph 700. Referring toFIG. 8B, the reference haplotype graph is extended to permit genotypingerrors and a new edge 852 is added to the graph, permitting a reading of“G” instead of “T” at this locus. The probability associated with thisedge is determined based on the rate of genotyping error for thegenotyping technology used. The unphased genotype data can therefore besplit into haplotypes “ACGCGC” and “GTGTGG”, the latter being a new pathbased on the extended reference haplotype graph. In some embodiments, toaccount for genotyping error, the out-of-sample phaser explicitly allowsgenotyping error with a constant probability of γ (which depends on theerror rate of the given technology, and is set to 0.01 in some cases)for each emitted edge label.

The example graphs shown include a small number of nodes and edges, andthus represent short sequences of genotype data. In practice, the beginstate node corresponds to the first locus on the chromosome and theaccepting state node the last locus on the chromosome, and the number ofedges in a path corresponds to the number of SNPs in a chromosome (L),which can be on the order of 50,000 in some embodiments. The thickestportion of the graph (i.e., a locus with the greatest number of possiblepaths), which depends at least in part on the DNA sequences ofindividuals used to construct the graph (K), can be on the order of5,000 in some embodiments. A large number of computations would beneeded (O(LK⁴) in the worst case) for a naive implementation of adynamic programming solution based on the Viterbi algorithm.

In some embodiments, the paths are pruned at each state of the graph tofurther improve performance. In other words, only likely paths are keptin the modified graph and unlikely paths are discarded. In someembodiments, after i markers (e.g., 3 markers), paths with probabilitiesbelow a certain threshold E (e.g., less than 0.0001%) are discarded. Forexample, a haplotype along a new path that accounts for bothrecombination and switching error would have very low probability ofbeing formed, and thus can be discarded. As another example, in the caseof unphased genotype data of “AG, CT, GG, CT, GG, CG,” (SEQ ID NO: 2) anew haplotype accounting for recombination can be forged by switchingpaths several times along the graph (additional edges would need to beadded but are not shown in the diagram). Given the low probabilityassociated with each switch, however, the formation of such a haplotypeis very unlikely and would be pruned from the resulting graph, while thepath that includes the genotyping error 825 has sufficiently highprobability, and is kept in the graph and used to thread the unphasedgenotype data into phased genotype data. By pruning unlikely paths fromthe modified graph, the dynamic programming-based phasing process isprevented from exploring very unlikely paths in the graph when threadinga new haplotype along it. The choice of ϵ determines the trade-offbetween the efficiency of the algorithm (in both time and space) and therisk of prematurely excluding the best Viterbi path. Computation savingsprovided by pruning can be significant. In some cases, phasing using anaive implementation can require 15 days per person while phasing withpruning only requires several minutes per person.

In some embodiments, the nodes and edges of the haplography can berepresented as follows:

struct Node {    int32_t id;    int32_t level;    Edge *outgoing[2]; };struct Edge {    int16_t id;    int8_t allele;    float weight;    Node*to; }

Even with a pruned haplotype graph, the number of nodes and edges can belarge and using the above data structures to represent the graph wouldrequire a vast amount of memory (on the order of several gigabytes insome cases). In some embodiments, the graph is represented in acompressed form, using segments. The term “segment” used herein refersto the data structure used to represent the graph in a compressed formand is different from the DNA segments used elsewhere in thespecification. Each segment corresponds to a contiguous set of edges inthe graph, with the following constraints: the end of the segment has upto 1 branch (0 branches are permitted), and no segment points to themiddle of another segment. In some embodiments, the data structure of asegment is represented as follows:

struct Segment {    int32_t timestamp;    int32_t index;    int32_tbegin;    int32_t end;    int32_t count[2];    Segment *edges[2]; }

FIG. 9 is a diagram illustrating an embodiment of a compressed haplotypegraph with segments. In this example, dashed shapes are used toillustrate the individual segments enclosed within. In some cases, acompressed graph associated with a chromosome can be represented usingseveral megabytes of memory, achieving memory reduction by a factor of1000 compared to the naïve implementation of nodes and edges.

Trio-Based Phasing

On a system such as the personal genomics services platform provided by23andMe®, DNA sequence information of one or both parents of theindividual is sometimes available and can be used to further refinephasing. With the exception of sites where all three individuals areheterozygous, the parental origin of each allele can be determinedunambiguously. For ambiguous sites, knowledge of patterns of locallinkage disequilibrium can be used to statistically estimate the mostlikely phase. In some embodiments, a refinement process that accountsfor parental DNA sequence information, referred to as trio-basedphasing, is optionally performed following the population-based phasingprocess to correct any errors in the output of the population-basedphasing process and improve phasing accuracy. In some embodiments, thetrio-based phasing technique is a post-processing step to be applied tosequences for which a previous population-based linkage-disequilibriumphasing approach has already been applied. The trio-based phasingtechnique can be used in combination with any existing phasing processto improve phasing quality, provided that an estimate of the switcherror rate (also referred to as the phasing error rate) is available.

In some embodiments, trio-based phasing receives as inputs a set ofpreliminary phased haplotype data (e.g., output of an out-of-samplepopulation-based phasing technique described above), and employs aprobabilistic graphic model (also referred to as a dynamic Bayesiannetwork) that models the observed alleles, hidden states, andrelationships of the parental and child haplotypes. The input includesthe set of preliminary phased haplotype data as well as the phasedhaplotype data of at least one parent. The genotype data at a particularsite (e.g., the i-th SNP on a chromosome) for each individual in thetrio (i.e., mom, dad, or child (i.e. the individual whose genetic datais being phased)) are represented by the following variables:

G*₀ ^(,i), G*₁ ^(,i) ∈{0, 1}: the observed alleles for haplotypes 0 and1, provided as input data. For the child, the input data can be obtainedfrom the output of the population-based phasing process (e.g., thepreliminary haplotype data). For the parent, the input data can be theoutput of the population-based phasing process or the final output of arefined process.

H*_(m) ^(,i), H*_(p) ^(,i) ∈{0, 1}: the hidden true alleles of theindividual's maternal (m) and paternal (p) haplotypes.

P*^(,i) ∈{m, p}: a hidden binary phase indicator variable that is set tom whenever G*₀ ^(,i) corresponds to H*_(m) ^(,i) and set to p wheneverG*₀ ^(,i) corresponds to H*_(p) ^(,i).

The relationship between parental and child haplotypes are encoded bytwo additional variables, T^(mom,i), T^(dad,i) ∈{a, b}, where aindicates transmission of the parent's maternal haplotype to the childand b indicates transmission of the parent's paternal haplotype to thechild. In some embodiments, a=0 and b=1.

The following assumptions are made about the model:

-   -   1. The hidden true alleles for each parent at each position        (i.e., H_(*) ^((mom,dad),i)), the initial phase for each        individual (i.e., P*^(,1)), and the initial transmission for        each parent (i.e., T*^(,1)) are independently drawn from uniform        Bernoulli priors.    -   2. The phase indicator variables for each individual and the        transmission indicator variables for each parent are each        sampled according to independent first order Markov processes.        Specifically,

${P\left( {P^{*{,i}}❘P^{*{,{i - 1}}}} \right)} = \left\{ \begin{matrix}{{1 - {s{if}P^{*{,i}}}} = P^{*{,{i - 1}}}} \\{s{otherwise}}\end{matrix} \right.$${P\left( {T^{*{,i}}❘T^{*{,{i - 1}}}} \right)} = \left\{ \begin{matrix}{{1 - {r{if}T^{*{,i}}}} = T^{*{,{i - 1}}}} \\{r{otherwise}}\end{matrix} \right.$

where s is the estimated switch error probability between consecutivesites in the input haplotypes and r is the estimated recombinationprobability between sites in a single meiosis. In some embodiments, s isset to a default value of 0.02 and r is set to a default value of

${\frac{1}{2}\left( {1 - e^{{- 2}{(\frac{37.5}{500000})}}} \right)} \approx {0.000075.}$

-   -   3. The hidden true alleles for the child at each position (i.e.,        H_(*) ^(kid,i)) are deterministically set on the parents' true        hidden haplotypes (i.e., neglecting the possibility of private        mutations) and their respective transmission variables.    -   4. The observed alleles are sampled conditionally on the true        alleles and the phase

variables with genotyping error, according to the following model:

${P\left( {{G_{0}^{*{,i}}❘H_{m}^{*{,i}}},H_{p}^{*{,i}},P^{*{,i}}} \right)} = \left\{ \begin{matrix}{{1 - {g{if}G_{0}^{*{,i}}}} = H_{p^{*{,t}}}^{*{,i}}} \\{g{otherwise}}\end{matrix} \right.$

according to the estimated genotyping error rate.

The following expression is used to characterize the trio-based phasingprocess:

maximize Pr(H _(m) ^(kid) ,H _(p) ^(kid) ,H _(m) ^(mom) ,H _(p) ^(mom),H _(m) ^(dad) , H _(p) ^(dad))

given H ^(i) _(m) +H ^(i) _(p) =G _(p) ^(i) +G ₁ ^(i) ∀i∈{kid,mom,dad}

FIG. 10 is a diagram illustrating an embodiment of a dynamic Bayesiannetwork used to implement trio-based phasing. The diagram depicts thestructure of the dynamic Bayesian network using plate notation. Roundedrectangles (also referred to as plates) such as 1002 and 1004 are usedto denote repeated structures in the graph model. Each plate correspondsto a position (e.g., the i-th marker) on the individual's chromosome. Inplate 1002 which corresponds to position i−1, variables which are notconnected to any variables from other plates (e.g., H_(m) ^(kid,i−1))are omitted from the diagram. Plate 1004 shows a detailed template forposition i ∈{1, 2, . . . , L}. As shown, nodes represent randomvariables in the model, and edges represent conditional dependencies.Shaded nodes (e.g., node 1006) represent random variables which areobserved at testing time, and nodes with thickened edges (e.g., node1008) represent variables which have dependencies across plates.

Trio-based phasing includes using the probabilistic model to estimatethe most probable setting of all unobserved variables, conditioned onthe observed alleles. In some embodiments, the most probable H variablesare determined using a standard dynamic programming-based technique(e.g., Viterbi). One can visualize the model as plates corresponding toi ∈{1, 2, . . . , L} being stacked in sequential order, and the pathsare formed by the interconnections of nodes on the same plate, as wellas nodes across plates.

FIG. 11 is a flowchart illustrating an embodiment of a process toperform trio-based phasing based on the model of FIG. 10 . Process 1100can be performed on a system such as 100 or 206, and can be used toimplement phasing module 302 to perform post-processing ofpopulation-based phasing. It is assumed that a model such as 1000 isalready established.

At 1102, emission probabilities are precomputed for each plate of model1000. In some embodiments, the emission probabilities, which correspondto the most likely setting for the H variables given the G, P, and Tvariables, are found using a dynamic programming (e.g., Viterbi) basedprocess. Referring to FIG. 10 , for a given position i, there are 2possible settings (0 or 1) for each for the variables P^(mom), P^(kid),P^(dad), T^(mom), T^(dad); there are two possible settings (0 or 1) foreach of the six H variables; and there are 3 possible settings (0, 1 ormissing) or each of the six G variables, 2⁵*2⁶*3⁶=1.5 million possiblecombinations. In subsequent steps, a dynamic programming process willsearch these combinations to identify the most likely setting for the Hvariables.

At 1104, transition probabilities are computed based at least in part onthe values of transition probabilities from the previous position.Referring to FIG. 10 , at a given position i, the values of transitionvariables T and P are dependent on the values of the T and P variablesfrom the previous position. There are 2 possible settings (0 or 1) foreach of the 5 P and T variables in the upper box 1002, and 2 possiblesettings (0 or 1) for each of the 5 P and T variables in the lower box.The possible combinations of the T and P values are therefore2⁵*2⁵=1024.

At 1106, based on the computed probabilities, the settings of transitionvariables T and P across the entire chromosome sequence (i.e., for i=1,. . . , L) are searched to determine the settings that would most likelyresult in the observed values. In some embodiments, the determination ismade using a dynamic programming technique such as Viterbi, and 2⁵*2⁵*Lstates are searched.

At 1108, the setting of H variables is looked up across the entiresequence to determine the settings that would most likely result in thegiven G, P, and T variables. This requires L table lookups.

The trio-based phasing solves the most likely settings for the Hvariables (the hidden true alleles for the individual's maternal andpaternal haplotypes at a given location). The solution is useful forphasing the child's DNA sequence information as well as for phasing aparent's DNA sequence information (if the parent's DNA sequenceinformation is unphased initially). In the event that only one parent'sDNA sequence information is available, the other parent's DNA sequenceinformation can be partially determined based on the DNA sequenceinformation of the known parent and the child (e.g., if the child'salleles at a particular location is “AC” and the mother's alleles at thesame location are “CC”, then one of the father's alleles would be “A”and the other one is unknown). The partial information can be marked(e.g., represented using a special notation) and input to the model. Thequality of trio-based phasing based on only one parent's information isstill higher than population-based phasing without using the trio-basedmethod.

In addition to improved haplotypes data, the result of trio-basedphasing also indicates whether a specific allele is deemed to beinherited from the mother or the father. This information is stored andcan be presented to the user in some embodiments.

Local Classification

Local classification refers to the classification of DNA segments asoriginating from an ancestry associated with a specific geographicalregion (e.g., Eastern Asia, Scandinavia, etc.) or ethnicity (e.g.,Ashkenazi Jew).

Local classification is based on the premise that, T generations ago,all the ancestors of an individual were unadmixed (i.e., originatingfrom the same geographical region). Starting at T generation, ancestorsfrom different geographical regions produced admixed offspring. Geneticrecombination breaks chromosomes and recombines them at each generation.After T generations, 2T meiosis occurred. As a result, the expectedlength of a recombination-free segment is expressed as:

$L = {\frac{F}{2T}{cM}}$

where F corresponds to a fixed length segment. In some embodiments, theexpected length L is determined to be 100 SNPs. It is used as thesegment size (also referred to as the window size) used in localclassification.

FIG. 12 is a flowchart illustrating an embodiment of a localclassification process. Process 1200 can be performed on a platform suchas 200 or a system such as 300.

Initially, at 1202, a set of K local ancestries is obtained. In someembodiments, the specification of the local ancestries depends on theancestries of unadmixed individuals whose DNA sequence information isused as reference data. For example, the set of local ancestries can bepre-specified to include the following: African, Native American,Ashkenazi, Eastern Asian, Southern Asian, Balkan, Eastern European,Western European, Middle Eastern, British Isles, Scandinavian, Finnish,Oceanian, Iberian, Greek, Sardinian, Italian, and Arabic. Many otherspecifications are possible; for example, in some embodiments the set oflocal ancestries correspond to individual countries such as the UK,Ireland, France, Germany, Finland, China, India, etc.

At 1204, a classifier is trained using reference data. In this example,the reference data includes DNA sequence information of unadmixedindividuals, such as individuals who are self-identified or identifiedby the system as having four grandparents of the same ancestry (i.e.,from the same region), DNA sequence information obtained from publicdatabases such as 1KG, HGDP-CEPH, HapMap, etc. The DNA sequenceinformation and their corresponding ancestry origins are input into theclassifier, which learns the corresponding relationships between the DNAsequence information (e.g., DNA sequence segments) and the correspondingancestry origins. In some embodiments, the classifier is implementedusing a known machine learning technique such as a support vectormachine (SVM), a neural network, etc. A SVM-based implementation isdiscussed below for purposes of illustration.

At 1206, phased DNA sequence information of a chromosome of theindividual is divided into segments. In some embodiments, phased data isobtained using the improved phasing technique described above. Phaseddata can also be obtained using other phasing techniques such as BEAGLE.The length of the segments can be a predetermined fixed value, such as100 SNPs. It is assumed that each segment corresponds to a singleancestry.

At 1208, the DNA sequence segments are input into the trained classifierto obtain corresponding predicted ancestries. In some embodiments, theclassifier determines probabilities associated with the set of localancestries (i.e., how likely a segment is from a particular localancestry), and the ancestry associated with the highest probability isselected as the predicted ancestry for a particular segment.

In some embodiments, one or more SVMs are used to implement theclassifier. An SVM is a known type of non-probabilistic binaryclassifier. It constructs a hyper plane that maximizes the distance tothe closest training data point of each class (in this case, a classcorresponds to a specific ancestry). A SVM can be expressed using thefollowing general expression:

$\left\{ \begin{matrix}{{\min\limits_{w,\xi,b}\frac{1}{2}{w}^{2}} + {C{\sum\limits_{i}\xi_{i}}}} \\{{y_{i}\left( {{w*x_{i}} - b} \right)} \geq {1 - {\xi_{i}{\forall i}}}} \\{\xi_{i} \geq {0{\forall i}}}\end{matrix} \right.$

where w is the normal vector to the hyper plane, C is a penalty term(fixed), the ξ are slack variables, x_(i) represents the features of thedata point i to be classified, and y_(i) is the class of data point i.

FIG. 13 is a diagram illustrating how a set of reference data points isclassified into two classes by a binary SVM.

Since an SVM is a binary classifier and there are K (e.g., 18) classesof local ancestries to be classified, the classification can bedecomposed into a set of binary problems (e.g., should the sequences beclassified as African or Native American, African or Ashkenazi, NativeAmerican or Ashkenazi, etc.). One approach is the “one vs. one”technique where a total of

$\begin{pmatrix}K \\2\end{pmatrix}$

classifiers are trained and combined to form a single local ancestryclassifier. Specifically, there is one classifier configured todetermine the likelihood that a sequence is African or Native American,another to determine African or Ashkenazi, another to determine NativeAmerican or Ashkenazi, etc. During the training process, reference dataof DNA sequences and their corresponding ancestries is fed to the SVMfor machine learning. When an ancestry prediction for a DNA sequencesegment is to be made, each trained SVM makes a determination aboutwhich one of the ancestry pair the DNA sequence segment more likelycorresponds to, and the results are combined to determine which ancestryis most likely. Specifically, the ancestry that wins the highest numberof determinations is chosen as the predicted ancestry. Another approachis the “one vs. all” technique where K classifiers are trained.

Several refinements can be made to improve the SVM. For example, thenumber of unadmixed reference individuals can vary greatly per ancestralorigin. If 700 samples are from Western Europe but only 200 samples arefrom South Asia, the imbalance in the number of samples can cause theWestern European-South Asian SVM to “favor” the larger class. Thus, thelarger class is penalized to compensate for the imbalance according tothe following:

$\left\{ \begin{matrix}{{\min\limits_{w,\xi,b}\frac{1}{2}{w}^{2}} + {\sum\limits_{G}{C_{G}{\sum\limits_{i}\xi_{i}}}}} \\{{y_{i}\left( {{w*x_{i}} - b} \right)} \geq {1 - {\xi_{i}{\forall i}}}} \\{\xi_{i} \geq {0{\forall i}}}\end{matrix} \right.$ $C_{G} \propto \frac{1}{❘G❘}$

where w is the normal vector to the hyper plane, C_(G) is a penalty termfor class G, the ξ are slack variables, x_(i) represents the features ofthe data point i to be classified, and y_(i) is the class of data pointi.

Another refinement is to encode strings of SNPs according to thepresence or absence of features. One approach is to encode one featureat each SNP according to the presence or absence of the minor allele.Another approach is to take substrings of length 2 which have 4 featuresper position and which can be encoded based on their presence or absenceas 00, 01, 10, and 11. A more general approach is to use a window oflength L, and encode (L−k+1)·2^(k) features of length k according to thepresence or absence of the features.

The general approach is not always feasible for practicalimplementation, given that there are

$\sum\limits_{k = 1}^{100}{\left( {L - k + 1} \right) \times 2^{k}}$

features in a window of length L. With L=100, this number isapproximately 10³⁰, too large for most memory systems. Thus, in someembodiments, a modified kernel is used. In some embodiments, aspecialized string kernel is used that computes the similarity betweenany two given windows as the total number of substrings they share. Thisapproach takes into account that even very similar windows contain sitesthat have mutated, resulting in common subsequences along with deleted,inserted, or replaced symbols. Therefore, the specialized string kernelis a more relevant way of comparing the similarity between two 100 SNPwindows, and achieves much higher accuracy than the standard linearkernel.

Another refinement is to use supervised learning. Supervised learningrefers to the task of training (or learning) a classifier using apre-labeled data, also referred to as the training set. Specifically, anSVM classifier is trained (or learned) using a training set of customerswhose ancestry was known (e.g., self-reported ancestries). Parameters ofthe SVM classifier are adjusted during the process. The trainedclassifier is then used to predict a label (ancestry) for any newunlabeled data.

Error Correction

The results of the local classifier can contain errors. FIG. 14 is adiagram illustrating possible errors in an example classificationresult. In this example, a portion of an individual's DNA sequence isphased and locally classified. The classification result 1400 is shown,in which the top portion corresponds to the haplotype inherited from oneparent and the bottom portion corresponds to the haplotype inheritedfrom the other parent. Each segment of the individual's DNA is labeledaccording to the local classifier output as Finland, Ashkenazi, orBritish Isles. Two types of possible errors are illustrated. Starting atlocation 1402, the label switched from Finland to Ashkenazi for thefirst half of the chromosome and from Ashkenazi to Finland for thesecond half. This switch can be caused by a double recombination, orsimply be a phasing error. At locations 1404 and 1406, there are twopredictions for segments that correspond to British Isles ancestry.These repeated predictions can be evidence of ancestry changes, or becaused by correlated local classifier (e.g., SVM) prediction errors.Thus, error correction is implemented in some embodiments. In addition,because the output of the SVM tends to be noisy, the error correctionprocess also reconciles these local ancestry estimates and producescleaner classifications (i.e., smooths out the noise). In someembodiments, confidence scores for the ancestry assignments are alsocomputed.

In some embodiments, error correction is implemented using a HiddenMarkov Model (HMM), which is a statistical model in which the input datais being modeled as a Markov process with unobserved (hidden) states. Inan HMM, the observed signal (the input data) is being generated by ahidden process in a sequential manner. A standard HMI assumes that anobservation, given the hidden state that generated it, is independent ofall previous observations. The hidden state at any given position onlydepends on the hidden state at the previous position. In someembodiments, the input (observed data) to the HMI includes the predictedancestries of DNA sequence segments (e.g., the ancestries as predictedby the local classifier for segments that are 100 SNPs in length). Thehidden state corresponds to the true ancestries of the segments. Theoutput of the HMI forms a set of smoothed ancestry origins for thesegments. FIG. 15 is a graph illustrating embodiments of Hidden MarkovModels used to model phasing errors. In the example shown, the observedvariables O¹, O², O³, and O⁴ correspond to Finland, Finland, BritishIsles, and British Isles, respectively. The possible values for eachhidden state H correspond to the set of local ancestries. Given theobserved variables and the probabilities associated with the model, themost likely sequence of hidden states can be solved.

FIG. 16 shows two example graphs corresponding to basic HMMs used tomodel phasing errors of two example haplotypes of a chromosome. In thisexample, graph 1602 illustrates an HMI corresponding to one haplotypeand graph 1604 illustrates an HMI corresponding to another haplotype.The most likely sequences of the H variables have been solved and the Hvariables are labeled.

The basic model averages out the output of the learning machine andgenerates a smoother and less noisy output, but does not correct many ofthe errors in the output. For example, in FIG. 16 , two HMIs model twohaplotypes separately and independently. In practice, however, theswitching errors in two haplotypes are not independent but instead occurtogether. Thus, in some embodiments, the basic HMI is improved bytreating the true ancestries associated with two haplotypes at a givenlocation as a single hidden state. FIG. 17 shows an example graph inwhich one HMI jointly models both haplotypes to account for potentialphasing errors. This model is referred to as a Pair Hidden Markov Model(PHMM). In this case, the pair of true ancestries of both haplotypes ata given location is treated as the hidden state at the location. ThePHMM accounts for the phasing errors (e.g., error at location 1402 ofFIG. 14 ) that are present in the output of the local classifier.

Also, in FIG. 16 , the basic HMMs assume that the observed data arecomputed independently (in other words, the Os are uncorrelated). Inpractice, the classification errors are correlated. For example, sincelong stretches of DNA can be inherited from one generation to the next,the true ancestries of adjacent DNA segments can be correlated. Thus, insome embodiments, the model of FIG. 17 is further improved to model theinterdependencies between observed states. The improved model is anothertype of PHMM, referred to as an Autoregressive Pair Hidden Markov Model(APHMM). FIG. 18 shows an example graph of an Autoregressive Pair HiddenMarkov Model (APHMM). In this example, a given observed state O^(i)depends both on its corresponding underlying hidden state H^(i) and onthe previous observed state O^(i−1).

The graph defines a probabilistic model as follows:

Pr(H ¹ ,H ² ,H ³ , . . . , O ¹ ,O ² ,O ³, . . . ) =Pr(H ¹)Pr(O ¹ |H¹)Pr(H ² |H ¹)Pr(O ² |H ² ,O ¹)Pr(H ³ |H ²)Pr(O ³ |H ³ ,O ²) . . .

where probabilities P(O^(j)|H^(j), O^(j−1)) are referred to as theemission parameters, and probabilities Pr(H^(i)|H^(i−1)) are referred toas the transition parameters.

The model outputs probabilities associated with ancestry assignments ofthe most probable sequence. Training is required to estimate theemission parameters and the transition parameters. In some embodiments,an expectation maximization method is used to estimate the parameters.

The emission parameters characterize how well the local classifierpredicts the ancestry. Specifically, given the underlying true state ofa segment, what is the probability that the local classifier will outputthe true state. FIG. 19 is an example data table displaying the emissionparameters. The data table can be generated based on empirical data ofthe local classifier's output and reference data of unadmixedindividuals. For example, suppose there are 1000 individuals in thereference population who have identified themselves as unadmixedAshkenazis. The local classifier output, however, labels 1% of theseindividuals' segments as African, 3% as Native American, 40% asAshkenazi, etc. Thus, P (O=African|H=Ashkenazi) is 1%, P (O=NativeAmerican|H=Ashkenazi) is 3%, P (O=Ashkenazi|H=Ashkenazi) is 40%, etc.These values are used to fill the corresponding entries in the table.The values for emission parameters are found along the diagonal of thetable.

The transition parameters correspond to the probability of a particularhidden state of the model given the previous hidden state. Theyrepresent the statistical likelihood of observing certain sequences oftrue ancestries in the population, and therefore need to be determinedbased on admixed data. However, it is not possible to obtain fullytransitioned and accurately labeled genomes from actual admixedindividuals. Thus, to determine the transition parameters, an iterativeapproach is used. Initially, the transition parameters are arbitrarilychosen to establish an initial model. The initial model is used toperform error correction. Based on the error corrected results, themodel is updated by applying an expectation maximization method. Theprocess can be repeated until a final convergent model is achieved.

Once the emission parameters and the transition parameters areestablished, the model is fully specified. Thus, the most likelysequence of hidden variables can be determined based on the observedstates using conventional HMI techniques. For example, a probabilisticscoring scheme is used to determine the most likely sequence in someembodiments. All the possibilities associated with the hidden states arelisted, and a set of scoring rules are specified to reward or penalizecertain events by adding or subtracting a score associated with asequence. For example, a change in adjacent haplotypes is likely anerror; therefore, whenever two adjacent haplotypes are different, thescore is reduced according to the rules. A mismatched observedstate/hidden state pair also indicates likely error; therefore, wheneverthere is a mismatch of predicted ancestry and the underlying ancestry,the score is reduced. The most likely sequence of hidden states can bedetermined by computing scores associated with all possible combinationsof observed states and hidden states, and selecting the sequence thatleads to the highest score. In some embodiments, more efficienttechniques for selecting the most likely sequence such as dynamicprogramming are employed to break the problem down into subproblems andreduce the amount of computation required. For example, the problem canbe reduced to recursively determine the best ancestry assignment foreverything to the left or the right of a particular position.

As described above, training is required to obtain parameters for thePHMM (or APHMM). In some embodiments, an ensemble technique is usedwhere reference population is grouped into distinct subsets to serve asdifferent types of training data resulting in different types of models.For example, different types of reference individuals that tend to havesimilar ancestries are identified and grouped into subsets. Such subsetscan be formed from admixed individuals (e.g., Latinos, Eurasians, etc.).as well as unadmixed individuals (e.g., East Asians, Northern Europeans,etc.) Data from a subset is used to determine the parameters of themodel for that subset. The resulting model is a model specific to thereference group (e.g., a Latino-specific model, a Eurasian-specificmodel, an East Asian specific-model, a Northern European-specific model,etc.). In some embodiments, the error correction process applies itsinput to all available models, and the results are weighted based onconfidence measures and then combined according to a Bayesian modelaveraging technique.

FIGS. 20A-20D are ancestry assignment plots illustrating differentresults that are obtained using different techniques. FIG. 20A is anexample plot of ancestry assignments of two simulated haplotypes of achromosome of a simulated individual. Each pattern in the haplotyperepresents a corresponding ancestry origin. This plot is referred to asthe true ancestries of the haplotypes. This data is input into differenttypes of ancestry assignment modules and the outputs are compared.First, the data is input into a local classifier without errorcorrection to generate an output that corresponds to FIG. 20B. As can beseen, the output is very noisy, where an underlying region that actuallyonly corresponds to a single ancestry is given many ancestry assignmentsby the local classifier. The local classifier output is input into abasic HMM (e.g., the HMM of FIG. 16 ) that performs an averaging andcorrects for the noise in the local classifier output, one haplotype ata time, to generate an output that corresponds to FIG. 20C. The HMMcorrects some of the noise in the local classifier output, but does notcorrect phasing errors or correlated classifier errors, and thereforethe result does not perfectly resemble the true ancestry assignments ofthe two haplotypes. The local classifier output is also input into anHMM that corrects for noise, phasing errors and correlated classifiererrors (e.g., the PHMM of FIG. 18 or the APHMM of FIG. 19 ) to generatean output that corresponds to FIG. 20D. The output perfectly matches thetrue ancestries of the underlying haplotypes (even though the order ofthe haplotypes are reversed in the output).

FIG. 21 is a table comparing the predictive accuracies of ancestryassignments with and without error correction. The “before” columns showthe recall and the precision associated with unadmixed or admixedindividuals using a basic HMI. The “after” columns show the recall andthe precision associated with unadmixed or admixed individuals usingAPHMM. As used herein, recall refers to what proportion of a particularancestry is correctly predicted, and precision refers to what proportionof a particular ancestry prediction is correct. For example, if 20% ofthe underlying reference data corresponds to African ancestry, and theassignment technique predicts that a given haplotype segment correspondsto African ancestry 10% of the time, then recall refers to the portionof the 20% of the African ancestry that is correctly predicted, andprecision refers to the portion of the 10% of the African ancestrypredictions that in fact corresponds to true African ancestry. As can beseen, recall and prediction of both unadmixed individuals and admixedindividuals are improved post-error correction.

Recalibration

The error correction module outputs the most probable sequence ofancestry assignments for a pair of haplotypes, and posteriorprobabilities associated with the corresponding assignments. Theposterior probabilities are recalibrated to establish confidencemeasures associated with the ancestry assignments. A well calibratedprediction with a probability of P should be correct P of the times. Howwell the posterior probability of the output is calibrated can bedetermined based on reference data of actual unadmixed individualsand/or simulated admixed individuals. For example, if it is known thatin the reference data, 10% of the haplotype segments correspond to EastEuropean ancestry, but the output predicts with 80% posteriorprobability that 20% of all the haplotype segments correspond to EastEuropean ancestry, then the posterior probability is overly confident.By tallying the percentage of the reference data that corresponds to aspecific ancestry, and applying the reference data to the predictiveengine to obtain the posterior probability, a reliability plot ofaccuracy vs. posterior probability can be determined for each referencepopulation corresponding to a specific ancestry. FIG. 22A illustratesexample reliability plots for East Asian population and East Europeanpopulation before recalibration. In the example plots shown, lines 2202and 2204 indicate the ideal accuracy—posterior probabilitiescorrespondence, and lines 2206 and 2208 indicate the actualaccuracy—posterior probabilities correspondence. Other reliability plotscan be similarly determined.

In some embodiments, Platt's recalibration technique is used torecalibrate the posterior probabilities. Logistics regression is appliedto posterior probabilities. A feature matrix X (e.g., 2^(nd) degreepolynomials) is defined, and a fit is determined based on the following:

${\Pr\left( {y = {1❘X}} \right)} = \frac{1}{1 + {\exp\left( {X\theta} \right)}}$

K-class recalibration is then performed (K being the number of localancestries). In some embodiments, K logistic curves are fit andrenormalized. In some embodiments, K logistic curves are fit andmultinomial logistic regression (i.e., softmax) is performed accordingto the following:

${\Pr\left( {y = {i❘X}} \right)} = {\frac{\exp\left( {X\theta_{i}} \right)}{1 + {{\sum}_{k = 1}^{K - 1}{\exp\left( {X\theta_{i}} \right)}}}{\forall{i \leq K}}}$

FIG. 22B illustrates example reliability plots for East Asian populationand East European population after recalibration. It is worth noting inthe original East European plot of FIG. 22A, when the posteriorprobability exceeds 60%, the accuracy is very poor. Thus, in thecalibrated result for East European population shown in FIG. 22B,accuracy is clipped at 60%.

In some embodiments, an isotonic regression technique (e.g., theZadrozny and Elkan method) is used to recalibrate the posteriorprobabilities, where recalibrated probabilities are estimated aspercentages of well classified training examples falling in each bin.

Given the input of (y_(i), p_(i))_(i)=1, . . . , n, the input is sortedin increasing order of p_(i). ϕ_(i) that monotonically increases withp_(i) but close to y_(i) are found. In some embodiments, apool-adjacent-violators (PAV) algorithm is used to solve:

$\min\limits_{\phi_{1},...,\phi_{n}}{\sum\limits_{i = 1}^{n}\left( {y_{i} - \phi_{i}} \right)^{2}}$ϕ₁ ≤ ϕ₂ ≤ ⋯ ≤ ϕ_(n)

where y_(i) is the label predicted for individual i, p_(i) is theuncalibrated probability associated with the prediction and ϕ_(i) is therecalibrated probability.

In some embodiments, modified isotonic regression techniques are used.For example, p_(i) can be bracketed into bins, and weights proportionalto the bin sizes are introduced to reduce computational cost. As anotherexample, regularization terms can be introduced to ensure smoothness ofthe output curves as follows:

${\min\limits_{\phi_{1},...,\phi_{n}}{\sum\limits_{i = 1}^{n}{\omega_{i}\left( {y_{i} - \phi_{i}} \right)}^{2}}} + {C{\sum\limits_{i = 1}^{n - 1}\left( {\phi_{i} - \phi_{i + 1}} \right)^{2}}}$ϕ₁ ≤ ϕ₂ ≤ ⋯ ≤ ϕ_(n)

In some embodiments, separate calibration regimes are used forindividuals with different amounts of effective switch error.Specifically, separate calibration curves are fitted for unadmixedindividuals (who have a low rate of effective switching error) oradmixed individuals (who have a high rate of effective switching error).

Label Clustering

In some embodiments, the recalibrated results are required to meet athreshold level of confidence before they are presented to the user. Ifthe threshold level is unmet, the assignments are clustered and repeatedas necessary until a total confidence level meets the threshold level.FIG. 23A is a flowchart illustrating an embodiment of a label clusteringprocess. Process 2300 can be performed on a platform such as 200 or asystem such as 300. At 2302, calibrated posterior probabilities ofpredicted ancestries associated with a haplotype segment are obtained(for example, as output of the calibration module). At 2304, thecalibrated posterior probabilities are sorted according to their values.At 2306, it is determined whether the calibrated posterior probabilitywith the highest value meets the threshold. If so, the ancestryassociated with the highest calibrated posterior probability is deemedto be the ancestry of the haplotype segment and is output (e.g.,presented to the user) at 2308. If, however, the threshold is not met,then, at 2310, the probabilities are clustered (e.g., summed) to formone or more new probabilities associated with broader geographicalregions, and the new probabilities are clustered again at 2304.2304-2310 are repeated until the threshold is met and a predictedancestry of sufficiently high confidence is presented to the user.

FIG. 23B is an example diagram illustrating process 2300 of FIG. 23A. Inthis example, geographical locations associated with differentancestries are organized into a hierarchical tree. The root node of thetree is the world (a segment that is labeled the world means that thesegment is unclassified). The next level corresponds to the continents.The next level under the continents corresponds to subcontinents. Thenext level corresponds to individual countries or specific geographicalregions within each subcontinent. Other representations are possible.The calibrated posterior probabilities of predicted ancestriesassociated with the haplotype segment are as follows: 5% probability ofbeing Oceania, 60% probability of being British Isles, 20% probabilityof being West Europe, 15% probability of being Greece. Depending on thevalue set for the threshold, different predicted ancestries can bepresented. For example, if the probability is set at 60%, the ancestryof British Isles is deemed to be associated with the segment and ispresented to the user. If the threshold is set at 70%, then none of thecurrent ancestry labels meets the threshold level. Thus, theprobabilities associated with specific countries or regions areclustered into subcontinents and the probabilities are summed. As shown,the probabilities associated with the British Isles and West Europe areclustered to form a probability indicating that the segment is 80%likely to be North-West European in its origin. The segment wouldtherefore be labeled as North-West European. If, however, the thresholdis set at 90%, then even at this level, no probability associated with asingle node of the hierarchical tree meets the threshold. Thus, theprobabilities are clustered again, combining the probabilities ofNorth-West Europe and South Europe to form a probability that thesegment is 95% likely to be from Europe. Europe is then presented as thepredicted ancestry associated with the segment. As shown, if a segmentcannot meet the threshold at the continental level, it is deemed to beworld/generic.

The output of the label cluster outputs the predicted ancestry for eachhaplotype segment. In some embodiments, the information is stored in adatabase and/or sent to an application to be displayed.

Display of Ancestry Composition Information

In some embodiments, once the ancestries associated with theindividual's chromosomes are determined, the results can be presentedvia various user interfaces upon user requests. The user interfaces canalso present ancestry information obtained using other techniques solong as the data being presented includes requisite information such asthe specific ancestries and proportions of the individual's genotypedata that corresponds to the specific ancestries.

FIG. 24 is a flowchart illustrating an embodiment of a process fordisplaying ancestry information. Process 2400 can be performed on aplatform such as 200 or a device such as 100.

At 2402, a request to display ancestry composition of an individual isreceived. In some embodiments, the request is received from anapplication that allows users to interact with their genetic informationand/or ancestry data. Referring to FIG. 2 as an example, a user can makea request to display ancestry information via an application (e.g., astandalone application or a browser-based application) installed onclient device 202. The application provides user interfaces (e.g.,buttons, selection boxes, or other appropriate user interface widgets)as well as associated logic to interact with ancestry prediction system206 and/or process data. The user can be the individual or anotherperson with permission to view the individual's ancestry information.

Returning to FIG. 24 , at 2404, upon receiving the request and inresponse, ancestry composition information of the individual isobtained. The ancestry composition information includes informationpertaining to proportions of the individual's genotype data deemed tocorrespond to specific ancestries. In some embodiments, the ancestrycomposition information includes hierarchical information of differentgeographical regions (e.g., 40% of the individual's genome correspondsto European ancestry, of which 80% corresponds to North-Western Europeanancestry and 20% corresponds to Southern European ancestry; of theNorth-Western European ancestry, 30% corresponds to Finish ancestry and70% corresponds to Scandinavian ancestry; of the Southern Europeanancestry, 50% corresponds to Italian ancestry and 50% corresponds toGreek ancestry, etc.). In some embodiments, the ancestry compositioninformation is obtained directly from a source such as a database or theoutput of a pipelined ancestry prediction process. In some embodiments,raw data obtained from a source is further processed to obtain ancestrycomposition information. For example, if the raw data only includespredicted ancestry composition information per haplotype segment, thenthe segments and their ancestries are tallied to determine theproportions of the individual's genes that are deemed to correspond tothe specific ancestries (e.g., 1000 out of 5000 segments are deemed tocorrespond to Italian ancestry and therefore 20% of the individual'sgenes are deemed to correspond to Italian ancestry).

At 2406, the ancestry composition information is presented to bedisplayed via a user interface.

In some embodiments, the ancestry composition information is initiallydisplayed according to geographical regions and proportions ofancestries deemed to correspond to those geographical regions.Subsequently, the user can request different data to be displayed viauser interfaces provided by the application. A further user request isoptionally received at 2408. At 2409, the type of request is determined.In response to the further user request and depending on the type ofrequest, different information is displayed. As shown, if the userrequest is a request to display subregions of a specific ancestry,subregions and proportions of the individual's ancestries correspondingto the subregions are displayed (or caused to be displayed on a displaydevice by processors) at 2410 in response. If the user request is arequest to display ancestries inherited from one or more parents, suchinformation is displayed (or caused to be displayed) at 2412 ifavailable. If the user request is a request to display ancestrycomposition information for a specific chromosome, the proportions ofancestries associated with the specific chromosome is displayed (orcaused to be displayed) at 2414. Other types of requests/displays arepossible.

FIG. 25 is a diagram illustrating an embodiment of a regional view ofancestry composition information for an individual. In this example,proportions of the individual's autosomal chromosome segments thatcorrespond to specific continents are displayed. The information ispresented to the user using two different views: a circle view 2502displaying the proportions as sections on a circle, where differentvisual formats (e.g., colors and/or patterns) represent differentancestral continents; and a list view 2504 displaying the proportionsand corresponding continent names. Other views are possible (e.g., a mapview displaying the regions associated with the ancestries).

The user is provided with the ability to expand the regions and viewmore detailed information pertaining to subregions. FIG. 26 . is adiagram illustrating an embodiment of an expanded view of ancestrycomposition information for an individual. In this example, the user canrequest to expand the circle view by moving a cursor or other pointingmechanism over a section of the circle view (e.g., European section2602), or by clicking on an entry in the list view (e.g., European entry2604). In response, the application expands the section and/or the entryin the list to show subregions of ancestries for the individual. In someembodiments, the subregions are determined according to the hierarchicalinformation of the ancestry composition information. In the exampleshown, a new section 2604 is created to include three subregions ofEurope (Northwest Europe, South Europe, and Generic Europe) from whichthe individual's ancestries can be traced as well as proportions of theindividual's DNA attributed to these subregions, and additional entries2606 are created to list the subregions and the proportions.

The subregions can be further expanded. FIG. 27 is a diagramillustrating an embodiment of a further expanded view of ancestrycomposition information for an individual. In this example, the user canrequest to expand the circle view by moving the cursor or other pointingmechanism over a subregion section (e.g., Northwest Europe section2702), or by clicking on an entry in the list view (e.g., NorthwestEurope entry 2704). In response, the application expands the subregionsection and/or the entry to show more detailed composition. In thisexample, the Northwest European ancestry is shown to include ancestriesfrom British Isles, Scandinavia, and generic Northwest Europe. Theprocess can be repeated as long as more granular data is available.

In some embodiments, the ancestry composition information that isobtained includes inheritance information, including proportions of theindividual's ancestries that are deemed to be inherited from either thefather or the mother. In other words, the inheritance informationpertains to how much of the individual's DNA corresponding to a specificancestry is inherited from each parent. For example, the trio-basedphasing result can indicate that for chromosome 1, haplotype 0, segments1-20 correspond to Scandinavian ancestry and are inherited from themother, segments 21-45 correspond to Italian ancestry and are inheritedfrom the mother also, segments 46-73 correspond to Greek ancestry andare inherited from the father, and so on. The segments from either themother or the father and the corresponding ancestries of the segmentsare tallied, and proportions of the ancestries attributed to each parentare computed. The inheritance information computation can be donefollowing trio-based phasing, at the time the request to displayinheritance from parents is made, or at some other appropriate time.Ancestry composition information of how much of the individual's DNAcorresponding to a specific ancestry is inherited from each parent isdisplayed. FIG. 28 is a diagram illustrating an embodiment of aninheritance view. In this example, the view is divided into two areasshown side-by-side. Area 2802 shows ancestries inherited from the fatherand area 2804 shows ancestries inherited from the mother. Acontinent-level view of ancestry composition attributed to each parentis shown initially. For example, 49.1% of the individual's DNA is deemedto correspond to sub-Saharan African ancestry, of which 3.9% isinherited from the father and 45.2% is inherited from the mother. Theuser is provided with the options to selectively view subregions insimilar manners as described in connection with FIGS. 25-27 .

In some embodiments, since the ancestry deconvolution process is appliedto individual chromosomes and the results are stored on the basis ofindividual chromosomes, the user has the option to select a specificautosomal chromosome or an X-chromosome to view its ancestralcomposition. FIGS. 29-30 are diagrams illustrating embodiments of achromosome-specific view. In FIG. 29 , a list of the autosomalchromosomes (and/or X-chromosome) that has undergone ancestrydeconvolution is presented to the user. Only 8 chromosomes are displayedin the figure but more are available. The user has the option to selecta specific one to view its ancestral composition. Upon receiving theuser selection (e.g., chromosome 1), in FIG. 30 , ancestry compositionsassociated with the selected chromosome are displayed in response.Although a list view is shown, other views such as the circle view shownin the preceding diagrams can also be presented. Subregions andinheritance of ancestries from a particular parent can be shown insimilar manners as described in connection with FIGS. 25-28 .

A pipelined ancestry deconvolution process and display of results havebeen described. The accuracy of ancestry predictions is greatly improvedover existing techniques, and the results are presented in aninformative and user-friendly fashion.

Although the foregoing embodiments have been described in some detailfor purposes of clarity of understanding, the invention is not limitedto the details provided. There are many alternative ways of implementingthe invention. The disclosed embodiments are illustrative and notrestrictive.

1. A method, implemented using a computer comprising one or moreprocessors and memory, the method comprising: obtaining, from aclassifier and by the one or more processors, a plurality of initialgeographical ancestry classifications respectively associated withsegments of two phased haplotypes of a chromosome pair of an individual,each phased haplotype being inherited from one of two parents of theindividual, each segment of one of the two phased haplotypescorresponding to a segment of another of the two phased haplotypes;obtaining, by the one or more processors, a Pair Hidden Markov Model(PHMM) in which observed states correspond to ordered pairs of theinitial geographical ancestry classifications associated with twocorresponding segments of the two phased haplotypes, and in which hiddenstates correspond to ordered pairs of underlying geographical ancestryclassifications associated with the two corresponding segments of thetwo phased haplotypes; determining, by the one or more processors andusing the PHMM, two geographical ancestry classifications respectivelyassociated with the two phased haplotypes; correcting, by the one ormore processors, the initial geographical ancestry classifications usingthe two geographical ancestry classifications; and outputting, by theone or more processors, the initial geographical ancestryclassifications as corrected.
 2. The method of claim 1, whereindetermining the two geographical ancestry classifications comprises:determining a most likely sequence of the hidden states given theinitial geographical ancestry classifications.
 3. The method of claim 2,wherein the two geographical ancestry classifications include twosequences of geographical ancestry classifications, and wherein the mostlikely sequence of the hidden states comprises the two sequences ofgeographical ancestry classifications.
 4. The method of claim 2, whereindetermining the most likely sequence of hidden states comprisesdetecting a correlated prediction error that was caused by theclassifier.
 5. The method of claim 4, wherein the initial geographicalancestry classification as corrected rectifies the correlated predictionerror.
 6. The method of claim 1, wherein the PHMM is an AutoregressivePair Hidden Markov Model (APHMM) in which the observed states aredependent on their corresponding hidden states and previous observedstates.
 7. The method of claim 1, wherein correcting the initialgeographical ancestry classifications comprises: obtaining a pluralityof Pair Hidden Markov Models (PHMMs) including the PHMM, wherein each ofthe plurality of PHMMs corresponds to a distinct reference population;and determining the initial geographical ancestry classifications ascorrected based on a weighting of the plurality of PHMMs.
 8. The methodof claim 7, wherein determining the initial geographical ancestryclassifications as corrected comprises: applying the initialgeographical ancestry classifications associated with segments of thetwo phased haplotypes to the plurality of PHMMs, and performing aBayesian model averaging on outputs of the plurality of PHMMs.
 9. Themethod of claim 1, wherein obtaining the PHMM comprises: performing, bythe one or more processors on the initial geographical ancestryclassifications, dynamic programming to determine the PHMM.
 10. A systemcomprising one or more processors, and one or more memories configuredto store instructions that, when executed by the one or more processors,cause the system to perform operations comprising: obtaining, from aclassifier, a plurality of initial geographical ancestry classificationsrespectively associated with segments of two phased haplotypes of achromosome pair of an individual, each phased haplotype being inheritedfrom one of two parents of the individual, each segment of one of thetwo phased haplotypes corresponding to a segment of another of the twophased haplotypes; obtaining a Pair Hidden Markov Model (PHMM) in whichobserved states correspond to ordered pairs of the initial geographicalancestry classifications associated with two corresponding segments ofthe two phased haplotypes, and in which hidden states correspond toordered pairs of underlying geographical ancestry classificationsassociated with the two corresponding segments of the two phasedhaplotypes; determining, using the PHMM, two geographical ancestryclassifications respectively associated with the two phased haplotypes;correcting the initial geographical ancestry classifications using thetwo geographical ancestry classifications; and outputting the initialgeographical ancestry classifications as corrected.
 11. The system ofclaim 10, wherein determining the two geographical ancestryclassifications comprises: determining a most likely sequence of thehidden states given the initial geographical ancestry classifications.12. The system of claim 11, wherein the two geographical ancestryclassifications include two sequences of geographical ancestryclassifications, and wherein the most likely sequence of the hiddenstates comprises the two sequences of geographical ancestryclassifications.
 13. The system of claim 11, wherein determining themost likely sequence of hidden states comprises detecting a correlatedprediction error that was caused by the classifier.
 14. The system ofclaim 10, wherein correcting the initial geographical ancestryclassifications comprises: obtaining a plurality of Pair Hidden MarkovModels (PHMMs) including the PHMM, wherein each of the plurality ofPHMMs corresponds to a distinct reference population; and determiningthe initial geographical ancestry classifications as corrected based ona weighting of the plurality of PHMMs.
 15. The system of claim 14,wherein determining the initial geographical ancestry classifications ascorrected comprises: applying the initial geographical ancestryclassifications associated with segments of the two phased haplotypes tothe plurality of PHMMs, and performing a Bayesian model averaging onoutputs of the plurality of PHMMs.
 16. A non-transitorycomputer-readable medium storing program instructions that, whenexecuted by one or more processors of a computing system, cause thecomputing system to perform operations comprising: obtaining, from aclassifier, a plurality of initial geographical ancestry classificationsrespectively associated with segments of two phased haplotypes of achromosome pair of an individual, each phased haplotype being inheritedfrom one of two parents of the individual, each segment of one of thetwo phased haplotypes corresponding to a segment of another of the twophased haplotypes; obtaining a Pair Hidden Markov Model (PHMM) in whichobserved states correspond to ordered pairs of the initial geographicalancestry classifications associated with two corresponding segments ofthe two phased haplotypes, and in which hidden states correspond toordered pairs of underlying geographical ancestry classificationsassociated with the two corresponding segments of the two phasedhaplotypes; determining, using the PHMM, two geographical ancestryclassifications respectively associated with the two phased haplotypes;correcting the initial geographical ancestry classifications using thetwo geographical ancestry classifications; and outputting the initialgeographical ancestry classifications as corrected.
 17. Thenon-transitory computer-readable medium of claim 16, wherein determiningthe two geographical ancestry classifications comprises: determining amost likely sequence of the hidden states given the initial geographicalancestry classifications.
 18. The non-transitory computer-readablemedium of claim 17, wherein the two geographical ancestryclassifications include two sequences of geographical ancestryclassifications, and wherein the most likely sequence of the hiddenstates comprises the two sequences of geographical ancestryclassifications.
 19. The non-transitory computer-readable medium ofclaim 17, wherein determining the most likely sequence of hidden statescomprises detecting a correlated prediction error that was caused by theclassifier.
 20. The non-transitory computer-readable medium of claim 16,wherein correcting the initial geographical ancestry classificationscomprises: obtaining a plurality of Pair Hidden Markov Models (PHMMs)including the PHMM, wherein each of the plurality of PHMMs correspondsto a distinct reference population; and determining the initialgeographical ancestry classifications as corrected based on a weightingof the plurality of PHMMs.